package com.salinity.kun.algorithm.harmonic;

import java.text.DecimalFormat;

public class GaussianElimination {

	/**
	 * 将系数全部置于等式左边，右边为0 如{1,1,1,6} {1,3,-2,1} {2,-2,1,1},
	 * https://blog.csdn.net/cheneyshark/article/details/78751550 有两个限制：
	 * 每次运算时，必须保证对角线上的元素不为0（即运算中的分母不为0），否则算法无法继续进行。
	 * 即使a(kk)的值不为零，但如果绝对值很小，由于第k次运算中a(kk)在分母位置，因此作除数会引起很大的误差，从而影响算法的稳定性。
	 * 
	 * 
	 * @param matrix
	 * @return
	 */
	public static double[] sequentialElimination(double[][] matrix) {
		int rows = matrix.length;
		int cols = matrix[0].length;
		double[][] m = new double[rows][cols];
		int i;
		for (i = 0; i < rows; i++) {
			for (int j = 0; j < cols; j++) {
				m[i][j] = matrix[i][j];
			}
		}

		for (i = 0; i < m.length; i++) {
			doElimination(m, i);
		}
		return getRslt(m);
	}

	/**
	 * 列主消元
	 * 
	 * @param matrix
	 * @return
	 */
	public static double[] columnMainElimination(double[][] matrix) {

		int rows = matrix.length;
		int cols = matrix[0].length;
		double[][] m = new double[rows][cols];
		int i;
		// 拷贝数组
		for (i = 0; i < rows; i++) {
			for (int j = 0; j < cols; j++) {
				m[i][j] = matrix[i][j];
			}
		}

		int maxCol;
		for (i = 0; i < rows; i++) {

			// 从当前行开始查找
			maxCol = findMaxAbsRow(m, i, i);
			if (i != maxCol) {
				swapCol(m, i, maxCol);
			}
			doElimination(m, i);
		}
		return getRslt(m);
	}


	////////////////////////// private///////////////////

	private static void swapCol(double[][] m, int col1, int col2) {
		double temp;
		int n = m[0].length;
		for (int j = 0; j < n; j++) {
			temp = m[col1][j];
			m[col1][j] = m[col2][j];
			m[col2][j] = temp;
		}
	}

	/**
	 * 找到当前列最大绝对值最大的数所在的行
	 * 
	 * @param m
	 * @param column
	 * @param startRow
	 *            起始行
	 * @return
	 */
	private static int findMaxAbsRow(double[][] m, int column, int startRow) {
		int col = column;
		int r = m.length;
		int c = m[0].length;
		if (col >= c) {
			throw new RuntimeException("column is out of range.");
		}

		double maxVal = Integer.MIN_VALUE;
		int rsltCol = startRow, i = startRow;
		for (; i < r; i++) {
			if (maxVal < Math.abs(m[i][col])) {
				maxVal = Math.abs(m[i][col]);
				rsltCol = i;
			}
		}
		return rsltCol;
	}

	/**
	 * 回代
	 * 
	 * @param afterElimination
	 * @return
	 */
	private static double[] getRslt(double[][] afterElimination) {
		double[][] m = afterElimination;
		int r = m.length;
		int c = m[0].length;
		int i, j;
		double rslt[] = new double[r];
		// 回代,从倒数第二行开始,只处理最后一列
		rslt[r - 1] = m[r - 1][c - 1];
		for (j = r - 2; j >= 0; j--) {
			for (i = 0; i <= j; i++) {
				m[i][c - 1] -= m[i][j + 1] * m[j + 1][c - 1];
			}
			rslt[j] = m[j][c - 1];
		}
		return rslt;
	}

	/**
	 * 消元,处理当前行与余下的行
	 * 
	 * @param m
	 * @param ri
	 */
	private static void doElimination(double[][] m, int ri) {
		int i, j, k = ri, r = m.length, c = m[0].length;
		double temp;

		// 将主元行[i][i] 处系数化为1，即数据同除以[i][i]
		temp = m[k][k];
		for (j = k; j < c; j++) {
			m[k][j] /= temp;
		}

		// 将主元行分别乘以对应行的 系数（可以消掉当前列的未知数），再跟对应行求差
		for (i = k + 1; i < r; i++) {

			// 当前行所乘 的倍数，因为主元行系数已经化为1了，这里直接乘以当前行该列的系数
			temp = m[i][k];
			for (j = k; j < c; j++) {
				// 对应行的各列 - temp * 当前行的各列
				m[i][j] -= temp * m[k][j];
			}
		}
	}

	/////////////////////// main/////////////////////////
	public static void main(String[] args) {

		double[][] input = { { 1, 1, 1 }, { 1, 3, -2 }, { 2, -2, 1 } };
		double[] b = { 6, 1, 1 };
		double[][] input2 = { { 1, 1, 1, 6 }, { 1, 3, -2, 1 }, { 2, -2, 1, 1 } };
		// double[][] input2 = { { 2, 3, 11, 5,2 }, { 1, 1, 5, 2,1 }, {
		// 2,1,3,2,-3},{1,1,3,4,-3} };
		double[][] input3 = { { -0.002, 2, 2, 0.4 }, { 1, 0.78125, 0, 1.3816 }, { 3.996, 5.5625, 4, 7.4178 } };
		// printMatrix(input);
		// double[] rslt = sequentialElimination(input3);
		double[] rslt = columnMainElimination(input3);
		double val = 0.0;
		int i = 0;
		for (double d : rslt) {
			System.out.printf("%4f ", d);
			val += input3[2][i++] * d;
		}
		System.out.println("\n" + val);
		System.out.println();
		rslt = sequentialElimination(input3);
		val = 0;
		i = 0;
		for (double d : rslt) {
			System.out.printf("%4f ", d);
			val += input3[2][i++] * d;
		}
		System.out.println("\n" + val);
		// double X[] = sequentialElimination(input,b);
		// System.out.println("rslt: "+Arrays.toString(X));
		/*
		 * 书上高斯消元的例子： 1 1 1 1 3 -2 2 -2 1
		 * 
		 * 6 1 1
		 *
		 *
		 * 书上列主消元的例子： -0.002 2 2 1 0.78125 0 3.996 5.5625 4
		 * 
		 * 0.4 1.3816 7.4178
		 */
	}
}
